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Abstract. Following a recent work by Robnik and Romanovski (J.Phys.A: 
Math.Gen. 39 (2006) L35, Open Syst. & Infor. Dyn. 13 (2006) 197-222) we derive 
the explicit formula for the universal distribution function of the final energies in a 
^ \ time-dependent ID harmonic oscillator, whose functional form does not depend on the 

■ details of the frequency w(t), and is closely related to the conservation of the adiabatic 

invariant. The normalized distribution function is P(x) = ir~ 1 (2[i 2 — x 2 )~z, where 
OC ' x = Ex — E\ , E\ is the final energy, E\ is its average value, and /z 2 is the variance of 

Ei. Ei and /z 2 can be calculated exactly using the WKB approach to all orders. 

o 



PACS numbers: 05.45.-a, 45.20.-d, 45.30,+s, 47.52. +j 



In a recent work P3 12] Robnik and Romanovski have studied the energy evolution 
in a general ID time-dependent harmonic oscillator and studied the closely related 
questions of the conservation of adiabatic invariants EJ Ej • Starting with the 
ensemble of uniformly distributed (w.r.t. the canonical angle variable) initial conditions 
on the initial invariant torus of energy E , they have calculated the average final energy 
Ei, the variance /i 2 and all higher moments. The even moments are powers of /z 2 , whilst 
the odd moments are exactly zero, because the distribution function P{E\) of the final 
energies E% is an even function w.r.t. E\. In this Letter we derive explicit formula for 
P(Ei), namely we shall derive 

where x — E\ — Ei, and Re denotes the real part, so that (JTJ) is zero for \x\ > fiy/2. 
This we do by using the exact results for the higher (even) moments and by employing 
the characteristic function f(y) of P{x). 
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The dynamics of our system is described by the Newton equation 

q + cu 2 (t)q = (2) 

which is generated by the system's Hamilton function H = H(q,p,t), whose numerical 
value E(t) at time t is precisely the total energy of the system at time t, and in case of 
ID harmonic oscillator this is 

where q,p,M,u> are the coordinate, the momentum, the mass and the frequency of the 
linear oscillator, respectively. 

The dynamics is linear in q,p, as described by (j2J), but nonlinear as a function of 
Ld(t) and therefore is subject to the nonlinear dynamical analysis. By using the index 
and 1 we denote the initial (t = t ) and final (t = ti) value of the variables. 

The transition map $ maps initial conditions {qo,Po) onto the final conditions 

-(:M:H::)(:) < 4 > 

with det($) = ad — be = 1, and a,b,c,d can be calculated as shown in [TJ Let 
E = H(q ,p ,t = t ) be the initial energy and E\ = H(qi,pi,t = ti) be the final 
energy. Introducing the new coordinates, namely the action I = E/uj and the angle 0, 
and assuming the uniform distribution of initial angles <fi over the period 2tt, we can 
immediately calculate the final average energy E\ and the variance 
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E J ul 



(5) 



It is shown in reference j2] that for arbitrary positive integer m, we have (E\ — Ei) 2m ~ l 
and 

(2m- 1) 



{El _ El) 2m = y ^iu El _ Ex y\ . (6 ) 

ml V / 

Thus 2m-th moment of P{E\) is equal to (2m — l)!!// 2m /m!, and therefore, indeed, all 
moments of P{Ei) are uniquely determined by the first moment E%. Obviously, P{E\) 
is in this sense universal, because it depends only on the average final energy E\ and 
the ratio u>i/u of the final and initial frequencies, and does not depend otherwise on 
any details of tu(t). It has a finite support (E min , E max ) = (E\ — fi\/2,Ei + /iV 7 ^), it 
is an even distribution w.r.t. E\ = (E min + E max )/2, and has an integrable singularity 
of the type 1/a/x at both E min and E max . This singularity stems from a projection of 
the final ensemble at t± onto the curves of constant final energies E\ of H(q,p,ti). Of 
course, all that we say here for the distribution of energies E\ holds true also for the 
final action, the adiabatic invariant l\ = Ei/u\. It is perhaps worthwhile to mention 
that the moments of our distribution according to © grow as 2 m / ^7cm, whilst e.g. in 
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the Gaussian distribution they grow much faster, namely as 2 m T(m + 1/2)/^, where 
T(x) denotes the gamma function. 

Now we derive the distribution function (0) using the characteristic function f(y) 
of P(x), namely 



f(y) = / e*y*P(x)dx. (7) 

J — oo 

We see immediately that the n-th derivative at y = is equal to 

/oo 
{ix) n P{x)dx = i n a n , (8) 
-oo 

where a n is the n-th moment of P(x) (in particular o~ 2 = A* 2 ), i-e. 

✓■oo 

a n = x n P(x)dx. (9) 



Using the Taylor expansion for f(y) = Yl™=o ^-y^-y" and the expressions © we see at 



once 

00 :2m,. 2m 



m!(2m)! 

m=0 v 7 

and using the formula (2m — 1)!! = (2m)!/(2 m m!), we obtain, 

m=0 x / \ / 

which can be summed and is equal to the Bessel function jH] 

f(y) = J (wV2). (12) 

Knowing the characteristic function ()7ll2j) we now only have to invert the Fourier 
transform, namely 

P(x) = — / e-^f(y)dy = — / e^MV^^dy, (13) 

which is precisely equal to ((TJ) for \x\ < \i\[2 and is zero otherwise. (See reference |Hj.) 
Indeed, the distribution ((TJ) is normalized to unity as it must be. It is essentially the so- 
called (3(1/2, 1/2) distribution or also termed arc sine density after shifting the origin 
of x to 1/2 and rescaling of x. It is obvious that the distribution function does not depend 
on any details of u(t) and is in this sense universal for ID time-dependent harmonic 
oscillator in case of uniform canonical ensemble of initial conditions, i.e. uniform w.r.t. 
the canonical angle. 

The rigorous method of calculating E\ and /i 2 is explained in |IJ E] . In the general 
case when (j2J cannot be solved exactly, the WKB method jTUj can be applied to all 
orders, as has been explained in detail in j2]. 
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Let us now present the direct algebraic derivation of the energy distribution function 
(dJ. By definition we have 



i=i 



(14) 



where we have to sum up contributions from all four branches of the function <p(Ei). 
Let us denote x = E\ — E 1 , so that we have 

x = E (5cos(20) + 7sin(20)) = /a/2sui(20 + (15) 

where 5 and 7 are expressed in terms of a,b,c,d as shown in PQ 12], ^ ne var i ance 
is /i 2 = ^(tf 2 + 7 2 ), and t&nip = 5/7, so that <\> = \ arcsin — ^. Therefore 



I #2.1 



for all four solutions z = 1, 2, 3, 4 and from (JT3)l we get (JTJ) at once. 



The first derivation demonstrates the power of the approach employing the 
characteristic function f(y) of P(x), giving us new insights, whereas the second one 
leads elegantly to the final result, eliminating many parameters by a rather elementary 
substitution. In fact, this second method is geometrically obvious by the following 
argument: In phase space, the initial distribution lies on an ellipse. The final distribution 
lies on a different ellipse with the same area, with the points equidistributed in 
the canonical angle variable. This final ellipse intersects the energy contours of the 
Hamiltonian, and thus after squeezing the final Hamiltonian so that its energy contours 
are circles, the desired probability distribution is simply the distribution of radii of those 
circles through which the final ellipse passes. 

In nonlinear systems the entire theory expounded in PI |2] must be reformulated 
and as such it is an important open problem [TI]. For the case of a separatrix crossing 
some interesting numerical results have been obtained in ^2j, namely P{E\) there has 
a substantial structure and is by far not so simple as 
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